module datatype
    implicit none
    integer, parameter:: pre = 8
    integer::i, j, k, iter
    integer::m, n, nn, mm, maxstp, if_cfd, maxstp_cfd, mmm, nm
    integer::numprocs
    real, parameter:: pi = 3.141592653, v_inf = 100.0
    real(pre), parameter::c = 1.0, dkc = 1.0, dit = 1.0 !< cord line
    real(pre), allocatable::theta(:), error(:, :, :)
    real(pre), allocatable::alpha(:, :), beta(:, :), gamma(:, :), jacobian(:, :)
    real(pre):: r, max_res(3), start_time, end_time, t1, t2

    type cell_type
        real(pre):: kc, it
        real(pre), pointer:: x(:), y(:)
        real(pre), pointer::res(:)
        real(pre), pointer::phi(:)
        real(pre)::u, v, cp
    end type
    type(cell_type), pointer::cell(:, :)

end module

subroutine grid_procs
    use datatype
    use omp_lib
    implicit none

    write (*, '(A21)', advance='no') 'Input dimension n, m:'
    read (*, *) n, m !< n can't be much more bigger than m
    write (*, '(A25)', advance='no') 'Input max iteration step:'
    read (*, *) maxstp
    write (*, '(A33)', advance='no') ' "1" for CFD, "0" for only grid:'
    read (*, *) if_cfd
    if (if_cfd == 1) then
        write (*, '(A25)', advance='no') 'Input cfd iteration step:'
        read (*, *) maxstp_cfd
    end if
    write (*, '(A15)') 'Input CPU cores:'
    read (*, *) numprocs

    call cpu_time(start_time)
!$  start_time = omp_get_wtime()
    call bc_init
    call metric

    if (if_cfd == 1) then
        call solver
    end if

    return

end subroutine

subroutine bc_init
    use datatype
    implicit none
    !m = 60; n = 60;
    !maxstp = 10000
    !maxstp_cfd = 2
    !if_cfd = 1
    !numprocs = 3

    nm = m
    m = m + 1
    mmm = m
    r = 50*c
    allocate (cell(n, m)) !< points : n * m
    allocate (theta(m))
    do i = 1, n
        do j = 1, m
            allocate (cell(i, j)%x(1:2)); allocate (cell(i, j)%y(1:2)); allocate (cell(i, j)%res(1:3))
        end do
    end do
    theta = 0.0

    !< bc
    do i = 1, mmm
        !< wall
        theta(i) = 2.0*pi/nm*(i - 1) !< deg 2 rad
        cell(1, i)%x(1) = 0.5*c*(1 + cos(theta(i)))
        if (theta(i) <= pi) then
            cell(1, i)%y(1) = 0.75*(0.2843*cell(1, i)%x(1)**3 - 0.1015*cell(1, i)%x(1)**4 - &
                                    0.3576*cell(1, i)%x(1)**2 - 0.1221*cell(1, i)%x(1) + 0.2969*cell(1, i)%x(1)**0.5)
        else
            cell(1, i)%y(1) = 0.75*(0.2843*cell(1, i)%x(1)**3 - 0.1015*cell(1, i)%x(1)**4 - &
                                    0.3576*cell(1, i)%x(1)**2 - 0.1221*cell(1, i)%x(1) + 0.2969*cell(1, i)%x(1)**0.5)*(-1.0)
        end if

        !< far
        cell(n, i)%x(1) = 0.5*r*cos(theta(i))
        cell(n, i)%y(1) = 0.5*r*sin(theta(i))

        cell(n, i)%x(2) = cell(n, i)%x(1)
        cell(n, i)%y(2) = cell(n, i)%y(1)
        cell(1, i)%x(2) = cell(1, i)%x(1)
        cell(1, i)%y(2) = cell(1, i)%y(1)

    end do

    !< inner init
    do i = 2, n - 1
        do j = 1, m
            cell(i, j)%x = cell(1, j)%x + (cell(n, j)%x - cell(1, j)%x)/(n - 1)*(i - 1)
            cell(i, j)%y = cell(1, j)%y + (cell(n, j)%y - cell(1, j)%y)/(n - 1)*(i - 1)
        end do
    end do

    return
end subroutine

subroutine metric
    use datatype
    use omp_lib
    implicit none
    nn = n - 1
    mm = m - 1
    max_res(1) = 0.0; max_res(2) = 0.0
    allocate (alpha(1:n, 1:m), beta(1:n, 1:m), gamma(1:n, 1:m))
    allocate (error(n, m, 3), jacobian(1:n, 1:m))
    error = 0.0
    !$omp parallel num_threads(numprocs)
    do iter = 1, maxstp!  !< loop
        !< periodic bc
        do i = 1, m
            do j = 1, n
                if (cell(j, i)%y(1) == 0.0 .and. cell(i, j)%x(1) >= c) then !< define cd
                    cell(j, 1)%x(1) = 0.5*(cell(j, 2)%x(1) + cell(j, m - 1)%x(1))
                    cell(j, 1)%y(1) = 0.5*(cell(j, 2)%y(1) + cell(j, m - 1)%y(1))
                    cell(j, m)%x(1) = cell(j, 1)%x(1)
                    cell(j, m)%y(1) = cell(j, 1)%y(1)
                end if
            end do
        end do

        do j = 1, m!< special work for wall &  far cells
            if (j == 1 .or. j == m) then !<四个顶点降为一阶精度
                alpha(1, 1) = (cell(2, 1)%x(1) - cell(1, 1)%x(1))**2 + &
                              (cell(2, 1)%y(1) - cell(1, 1)%y(1))**2

                beta(1, 1) = (cell(1, 2)%x(1) - cell(1, 1)%x(1))*(cell(2, 1)%x(1) - cell(1, 1)%x(1)) + &
                             (cell(1, 2)%y(1) - cell(1, 1)%y(1))*(cell(2, 1)%y(1) - cell(1, 1)%y(1))

                gamma(1, 1) = (cell(1, 2)%x(1) - cell(1, 1)%x(1))**2 + &
                              (cell(1, 2)%y(1) - cell(1, 1)%y(1))**2
!< ------------------------------------------------------------------------------------------------------
                alpha(n, 1) = (cell(n, 1)%x(1) - cell(n - 1, 1)%x(1))**2 + &
                              (cell(n, 1)%y(1) - cell(n - 1, 1)%y(1))**2

                beta(n, 1) = (cell(n, 2)%x(1) - cell(n, 1)%x(1))*(cell(n, 1)%x(1) - cell(n - 1, 1)%x(1)) + &
                             (cell(n, 2)%y(1) - cell(n, 1)%y(1))*(cell(n, 1)%y(1) - cell(n - 1, 1)%y(1))

                gamma(n, 1) = (cell(n, 2)%x(1) - cell(n, 1)%x(1))**2 + &
                              (cell(n, 2)%y(1) - cell(n, 1)%y(1))**2
!< ------------------------------------------------------------------------------------------------------
                alpha(1, m) = (cell(2, m)%x(1) - cell(1, m)%x(1))**2 + &
                              (cell(2, m)%y(1) - cell(1, m)%y(1))**2

                beta(1, m) = (cell(1, m)%x(1) - cell(1, m - 1)%x(1))*(cell(2, m)%x(1) - cell(1, m)%x(1)) + &
                             (cell(1, m)%y(1) - cell(1, m - 1)%y(1))*(cell(2, m)%y(1) - cell(1, m)%y(1))

                gamma(1, m) = (cell(1, m)%x(1) - cell(1, m - 1)%x(1))**2 + &
                              (cell(1, m)%y(1) - cell(1, m - 1)%y(1))**2
!< ------------------------------------------------------------------------------------------------------
                alpha(n, m) = (cell(n, m)%x(1) - cell(n - 1, m)%x(1))**2 + &
                              (cell(n, m)%y(1) - cell(n - 1, m)%y(1))**2

                beta(n, m) = (cell(n, m)%x(1) - cell(n, m - 1)%x(1))*(cell(n, m)%x(1) - cell(n - 1, m)%x(1)) + &
                             (cell(n, m)%y(1) - cell(n, m - 1)%y(1))*(cell(n, m)%y(1) - cell(n - 1, m)%y(1))

                gamma(n, m) = (cell(n, m)%x(1) - cell(n, m - 1)%x(1))**2 + &
                              (cell(n, m)%y(1) - cell(n, m - 1)%y(1))**2
!< ------------------------------------------------------------------------------------------------------
            else

                alpha(1, j) = (cell(2, j)%x(1) - cell(1, j)%x(1))**2 + &
                              (cell(2, j)%y(1) - cell(1, j)%y(1))**2

                beta(1, j) = (cell(1, j + 1)%x(1) - cell(1, j - 1)%x(1))*(cell(2, j)%x(1) - cell(1, j)%x(1))/2 + &
                             (cell(1, j + 1)%y(1) - cell(1, j - 1)%y(1))*(cell(2, j)%y(1) - cell(1, j)%y(1))/2

                gamma(1, j) = ((cell(1, j + 1)%x(1) - cell(1, j - 1)%x(1))/2)**2 + &
                              ((cell(1, j + 1)%y(1) - cell(1, j - 1)%y(1))/2)**2

                alpha(n, j) = (cell(n, j)%x(1) - cell(n - 1, j)%x(1))**2 + &
                              (cell(n, j)%y(1) - cell(n - 1, j)%y(1))**2

                beta(n, j) = (cell(n, j + 1)%x(1) - cell(n, j - 1)%x(1))*(cell(n, j)%x(1) - cell(n - 1, j)%x(1))/2 + &
                             (cell(n, j + 1)%y(1) - cell(n, j - 1)%y(1))*(cell(n, j)%y(1) - cell(n - 1, j)%y(1))/2

                gamma(n, j) = ((cell(n, j + 1)%x(1) - cell(n, j - 1)%x(1))/2)**2 + &
                              ((cell(n, j + 1)%y(1) - cell(n, j - 1)%y(1))/2)**2

            end if
        end do

        !$omp do
        do i = 2, n - 1
            do j = 2, m - 1
                alpha(i, j) = ((cell(i + 1, j)%x(1) - cell(i - 1, j)%x(1))/2)**2 + &
                              ((cell(i + 1, j)%y(1) - cell(i - 1, j)%y(1))/2)**2

                beta(i, j) = (cell(i, j + 1)%x(1) - cell(i, j - 1)%x(1))*(cell(i + 1, j)%x(1) - cell(i - 1, j)%x(1))/4 + &
                             (cell(i, j + 1)%y(1) - cell(i, j - 1)%y(1))*(cell(i + 1, j)%y(1) - cell(i - 1, j)%y(1))/4

                gamma(i, j) = ((cell(i, j + 1)%x(1) - cell(i, j - 1)%x(1))/2)**2 + &
                              ((cell(i, j + 1)%y(1) - cell(i, j - 1)%y(1))/2)**2

                cell(i, j)%x(2) = 0.5*(alpha(i, j)*(cell(i, j + 1)%x(1) + cell(i, j - 1)%x(1)) + &
                                       gamma(i, j)*(cell(i + 1, j)%x(1) + cell(i - 1, j)%x(1)) - &
           0.5*beta(i, j)*(cell(i + 1, j + 1)%x(1) + cell(i - 1, j - 1)%x(1) - cell(i - 1, j + 1)%x(1) - cell(i + 1, j - 1)%x(1))) &
                                  /(alpha(i, j) + gamma(i, j))

                cell(i, j)%y(2) = 0.5*(alpha(i, j)*(cell(i, j + 1)%y(1) + cell(i, j - 1)%y(1)) + &
                                       gamma(i, j)*(cell(i + 1, j)%y(1) + cell(i - 1, j)%y(1)) - &
           0.5*beta(i, j)*(cell(i + 1, j + 1)%y(1) + cell(i - 1, j - 1)%y(1) - cell(i - 1, j + 1)%y(1) - cell(i + 1, j - 1)%y(1))) &
                                  /(alpha(i, j) + gamma(i, j))

                cell(i, j)%res(1) = abs(cell(i, j)%x(2) - cell(i, j)%x(1))
                cell(i, j)%res(2) = abs(cell(i, j)%y(2) - cell(i, j)%y(1))
                error(i, j, 1) = cell(i, j)%res(1)
                error(i, j, 2) = cell(i, j)%res(2)

                cell(i, j)%x(1) = cell(i, j)%x(2) !< update
                cell(i, j)%y(1) = cell(i, j)%y(2)
            end do
        end do
        !$omp end do

        max_res(1) = maxval(error(1:n, 1:m, 1))
        max_res(2) = maxval(error(1:n, 1:m, 2))
        call cpu_time(end_time)
!$      end_time = omp_get_wtime()

        !$omp barrier
        !$omp single
        if (if_cfd == 0) then
            if (mod(iter, 5000) == 0 .or. iter == 1) then
                write (*, '(A55)') 'Max error(x)     Max error(y)     Iters     Time'
            end if
            if (mod(iter, 500) == 0) then
                write (*, '(2E20.12, I6, F10.2)') max_res(1), max_res(2), iter, end_time - start_time
            end if
        end if
        !$omp end single

    end do

    !< tec out
    !$omp single
    open (unit=111, file='test.plt', form='formatted')
    write (111, *) 'variables = "x", "y"'
    write (111, *) 'ZONE DATAPACKING = POINT,', 'I = ', m, 'J = ', n  !< i, j num must be included

    do i = 1, n
        do j = 1, m
            write (111, *) cell(i, j)%x(2), cell(i, j)%y(2)
        end do
    end do

    close (111)
    if (if_cfd == 1) then
        write (*, '(A19)') 'Grid process done !'
        write (*, '(A55)') 'Max error(x)     Max error(y)     Iters     Time'
        write (*, '(2E20.12, I6, F10.2)') max_res(1), max_res(2), iter - 1, end_time - start_time
    end if
    !$omp end single
    !$omp end parallel

    return
end subroutine

subroutine cfd_procs
    use datatype
    implicit none

    call cfd_

    stop
end subroutine

subroutine cfd_
    use datatype
    use omp_lib
    implicit none
    write (*, *)
    write (*, '(A19)') 'CFD process begin !'
    do i = 1, n
        do j = 1, m
            allocate (cell(i, j)%phi(1:2))
            cell(i, j)%u = 0.0
            cell(i, j)%v = 0.0
            cell(i, j)%cp = 0.0
            cell(i, j)%phi(:) = 0.0
        end do
    end do

    do i = 1, n
        do j = 1, m
            cell(i, j)%phi(1:2) = v_inf*cell(i, j)%x(2) !< init
        end do
    end do

    open (unit=123, file='cp.plt', form='formatted')
    !$omp parallel num_threads(numprocs)

    do iter = 1, maxstp_cfd !< loop

        !< periodic bc
        do i = 1, m
            do j = 1, n
                if (cell(j, i)%y(1) == 0.0 .and. cell(i, j)%x(1) >= c) then !< define cd
                    cell(j, 1)%phi(1) = 0.5*(cell(j, 2)%phi(1) + cell(j, m - 1)%phi(1))
                    jacobian(j, 1) = 0.5*(jacobian(j, 2) + jacobian(j, m - 1))
                    cell(j, m)%phi(1) = cell(j, 1)%phi(1)
                end if
            end do
        end do

        do j = 2, m - 1 !< special work for wall cells
            cell(1, j)%phi(2) = cell(2, j)%phi(1) - 0.5*beta(1, j)/gamma(1, j)* &
                                (cell(1, j + 1)%phi(1) - cell(1, j - 1)%phi(1))

            jacobian(1, j) = (cell(1, j + 1)%x(2) - cell(1, j - 1)%x(2))*(cell(2, j)%y(2) - cell(1, j)%y(2))/2 - &
                             (cell(1, j + 1)%y(2) - cell(1, j - 1)%y(2))*(cell(2, j)%x(2) - cell(1, j)%x(2))/2
            jacobian(1, j) = 1/jacobian(1, j)

            cell(1, j)%u = jacobian(1, j)*((cell(1, j + 1)%phi(1) - cell(1, j - 1)%phi(1))* &
                                           (cell(2, j)%y(2) - cell(1, j)%y(2)) - &
                                           (cell(2, j)%phi(1) - cell(1, j)%phi(1))*(cell(1, j + 1)%y(2) - cell(1, j - 1)%y(2)))/2

            cell(1, j)%v = jacobian(1, j)*((cell(2, j)%phi(1) - cell(1, j)%phi(1))* &
                                           (cell(1, j + 1)%x(2) - cell(1, j - 1)%x(2)) - &
                                           (cell(1, j + 1)%phi(1) - cell(1, j - 1)%phi(1))*(cell(2, j)%x(2) - cell(1, j)%x(2)))/2

            cell(1, j)%cp = 1 - (cell(1, j)%u**2 + cell(1, j)%v**2)/v_inf**2
        end do

        cell(1, 1)%phi(2) = 1/(1 - (beta(1, 1)/gamma(1, 1)))*(cell(2, 1)%phi(1) - (beta(1, 1)/gamma(1, 1))* &
                                                              cell(1, 2)%phi(1))
        jacobian(1, 1) = (cell(1, 2)%x(2) - cell(1, 1)%x(2))*(cell(2, 1)%y(2) - cell(1, 1)%y(2)) - &
                         (cell(1, 2)%y(2) - cell(1, 1)%y(1))*(cell(2, 1)%x(2) - cell(1, 1)%x(2))
        jacobian(1, 1) = 1/jacobian(1, 1)
!< -------------!-----------------------------------------------------------------------------------------
        cell(1, m)%phi(2) = 1/(1 - (beta(1, m)/gamma(1, m)))*(cell(2, m)%phi(1) + (beta(1, m)/gamma(1, m))* &
                                                              cell(1, m - 1)%phi(1))

        jacobian(1, m) = (cell(1, m)%x(2) - cell(1, m - 1)%x(2))*(cell(2, m)%y(2) - cell(1, m)%y(2)) - &
                         (cell(1, m)%y(2) - cell(1, m - 1)%y(1))*(cell(2, m)%x(2) - cell(1, m)%x(2))
        jacobian(1, m) = 1/jacobian(1, m)

        !$omp do
        do i = 2, n - 1
            do j = 2, m - 1

                cell(i, j)%phi(2) = 0.5*(alpha(i, j)*(cell(i, j + 1)%phi(1) + cell(i, j - 1)%phi(1)) + &
                                         gamma(i, j)*(cell(i + 1, j)%phi(1) + cell(i - 1, j)%phi(1)) - &
                               0.5*beta(i, j)*(cell(i + 1, j + 1)%phi(1) + cell(i - 1, j - 1)%phi(1) - cell(i - 1, j + 1)%phi(1) - &
                                                         cell(i + 1, j - 1)%phi(1)))/(alpha(i, j) + gamma(i, j))

                cell(i, j)%res(3) = abs(cell(i, j)%phi(2) - cell(i, j)%phi(1))
                error(i, j, 3) = cell(i, j)%phi(2) - cell(i, j)%phi(1)

                jacobian(i, j) = (cell(i, j + 1)%x(2) - cell(i, j - 1)%x(2))*(cell(i + 1, j)%y(2) - cell(i - 1, j)%y(2))/4 - &
                                 (cell(i, j + 1)%y(2) - cell(i, j - 1)%y(2))*(cell(i + 1, j)%x(2) - cell(i - 1, j)%x(2))/4

                jacobian(i, j) = 1/jacobian(i, j)

                !cell(i, j)%u = (cell(i, j +1)%phi(2) - cell(i, j - 1)%phi(2))/(cell(i, j +1)%x(2) - cell(i, j - 1)%x(2))
                !cell(i, j)%v = (cell(i +1, j)%phi(2) - cell(i - 1, j)%phi(2))/(cell(i + 1, j)%y(2) - cell(i - 1, j)%y(2))

                cell(i, j)%u = jacobian(i, j)*((cell(i, j + 1)%phi(1) - cell(i, j - 1)%phi(1))* &
                                               (cell(i + 1, j)%y(2) - cell(i - 1, j)%y(2)) - &
                                      (cell(i + 1, j)%phi(1) - cell(i - 1, j)%phi(1))*(cell(i, j + 1)%y(2) - cell(i, j - 1)%y(2)))/4

                cell(i, j)%v = jacobian(i, j)*((cell(i + 1, j)%phi(1) - cell(i - 1, j)%phi(1))* &
                                               (cell(i, j + 1)%x(2) - cell(i, j - 1)%x(2)) - &
                                      (cell(i, j + 1)%phi(1) - cell(i, j - 1)%phi(1))*(cell(i + 1, j)%x(2) - cell(i - 1, j)%x(2)))/4

                cell(i, j)%cp = 1 - (cell(i, j)%u**2 + cell(i, j)%v**2)/v_inf**2

            end do
        end do
        !$omp end do
        do i = 1, n
            do j = 2, m
                cell(i, j)%phi(1) = cell(i, j)%phi(2) !< update
            end do
        end do
        max_res(3) = abs(maxval(error(2:nn, 2:mm, 3)))

    end do

    call cpu_time(end_time)
!$  end_time = omp_get_wtime()
    !$omp barrier
    !$omp single

    write (*, '(A40)') 'Max error(phi)     Iters     Time'
    write (*, '(E20.12, I11, F10.2)') max_res(3), iter - 1, end_time - start_time
    write (*, 130)

    !do i = 2, n - 1
    write (123, *) 'variables = "x", "-Cp"'
    do j = 2, m
        write (123, *) cell(1, j)%x(2), -1.0*cell(1, j)%cp
    end do
    !end do
130 format(//9x, '--------- All computation is finished!---------'//)
!$omp end single
!$omp end parallel
    return
end subroutine

program gird
    use datatype
    implicit none

    call grid_procs

    stop
end program
